This script uses the output of r_l_preparation.Rmd and
produes all the values and plots reported in the paper.
if( !all(file.exists("./data/web_raw_trials.csv", "./data/web_experiment_cleaned.csv", "./data/field_raw_trials.csv", "./data/field_experiment_cleaned.csv", "./2_r_l_preparation.html")) )
{
# Seems like the script 2_r_l_preparation.Rmd was not invoked, so let's compile it!
require(rmarkdown);
#rmarkdown::render("./2_r_l_preparation.Rmd");
xfun::Rscript_call(rmarkdown::render, list(input="./2_r_l_preparation.Rmd"));
if( !all(file.exists("./data/web_raw_trials.csv", "./data/web_experiment_cleaned.csv", "./data/field_raw_trials.csv", "./data/field_experiment_cleaned.csv", "./2_r_l_preparation.html")) )
{
# Seems this issue is quite serious?
stop("Failed to compile the '2_r_l_preparation.Rmd' Rmarkdonw script: please try to do manually!\n");
}
}
# Load packages:
library(tidyverse) # data processing
library(brms) # bayesian models
#library(cmdstanr) # install it with: install.packages("cmdstanr", repos = c("https://mc-stan.org/r-packages/", getOption("repos"))) followed by install_cmdstan()
library(ggdist) # for plotting
library(tidybayes)
# option for Bayesian regression models:
# use all available cores for parallel computing
options(mc.cores = parallel::detectCores())
## Set the script's path as working directory
#parentfolder = rstudioapi::getActiveDocumentContext()$path
#setwd(dirname(parentfolder))
#parentfolder <- getwd()
parentfolder <- "."; # assume current folder is the document folder
#models <- paste0(parentfolder, '/models/')
plots <- paste0(parentfolder, '/plots/')
data <- paste0(parentfolder, '/data/')
if( !dir.exists("./cached_results") ) dir.create("./cached_results", showWarnings=FALSE);
# Various auxiliary functions:
library(parallel);
library(lme4);
library(performance);
library(brms);
library(bayestestR);
library(ggplot2);
library(gridExtra);
library(ggpubr);
library(sjPlot);
library(knitr);
brms_ncores <- max(detectCores(all.tests=TRUE, logical=FALSE), 4, na.rm=TRUE); # try to use multiple cores, if present
# Log odds to probability (logistic regression intercept):
lo2p <- function(x){ o <- exp(x); return (o/(1+o));}
# Log odds to odds ratio change as percent (logistic regression slopes):
lo2or <- function(x){ o <- exp(x); if(x < 0){ return (-(1-o)) } else { return (o-1) };}
# Log odds to odds ratio change between level values (logistic regression slopes):
#lo2ps <- function(a,b){if(b > 0){ return (plogis(a+b)-plogis(a)) }else{ return (-(plogis(a)-plogis(a+b))) }}
lo2ps <- function(a,b){ plogis(a+b) - plogis(a) }
# Scientific notation using Markdown conventions (inspired from https://www.r-bloggers.com/2015/03/scientific-notation-for-rlatex/):
scinot <- function(xs, digits=2, pvalue=TRUE)
{
scinot1 <- function(x)
{
sign <- "";
if(x < 0)
{
sign <- "-";
x <- -x;
}
exponent <- floor(log10(x));
if(exponent && pvalue && exponent < -3)
{
xx <- round(x / 10^exponent, digits=digits);
e <- paste0("×10^", round(exponent,0), "^");
} else
{
xx <- round(x, digits=digits+1);
e <- "";
}
paste0(sign, xx, e);
}
vapply(xs, scinot1, character(1));
}
# Escape * in a string:
escstar <- function(s)
{
gsub("*", "\\*", s, fixed=TRUE);
}
# Figure and Table caption adapted from https://stackoverflow.com/questions/37116632/rmarkdown-html-number-figures:
outputFormat = opts_knit$get("rmarkdown.pandoc.to"); # determine the output format of the document
if( is.null(outputFormat) ) outputFormat = ""; # probably not run within knittr
capTabNo = 1; capFigNo = 1; # figure and table caption numbering, for HTML do it manually
#Function to add the Table Number
capTab = function(x)
{
if(outputFormat == 'html'){
x = paste0("**Table ",capTabNo,".** ",x,"")
capTabNo <<- capTabNo + 1
}; x
}
#Function to add the Figure Number
capFig = function(x, show_R_version=TRUE, show_package_versions=NULL, is_map=FALSE)
{
if(outputFormat == 'html')
{
x <- paste0("**Figure ",capFigNo,".** ",x,"");
if( show_R_version || (!is.null(show_package_versions) && length(show_package_versions) > 0) )
{
x <- paste0(x, " Figure generated using ");
if( show_R_version ) x <- paste0(x, stringr::str_replace(R.version.string, stringr::fixed("R "), "[`R`](https://www.r-project.org/) "));
if( !is.null(show_package_versions) && length(show_package_versions) > 0 )
{
x <- paste0(x, ifelse( show_R_version, " and ", " "));
x <- paste0(x, ifelse( length(show_package_versions) > 1, "packages ", "package "));
x <- paste0(x, paste0(vapply(show_package_versions, function(x) paste0("`",x,"`", " (version ", packageVersion(x),")"), character(1)), collapse=", "), ".");
}
if( is_map ) x <- paste0(x, " Maps are using public domain data from the [Natural Earth project](https://www.naturalearthdata.com/) as provided by the `R` package `maps`.");
}
capFigNo <<- capFigNo + 1;
};
x;
}
# For reproducible reporting (see also the endof this document):
packageVersion('tidyverse')
packageVersion('ggplot2')
packageVersion('brms')
#packageVersion('cmdstanr')
packageVersion('ggdist')
# Load ggplot2 theme and colors:
source('theme_timo.R')
colorBlindBlack8 <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
# Load data:
web <- read_csv(paste0(data, 'web_experiment_cleaned.csv'))
web_raw <- read_csv(paste0(data, '/web_raw_trials.csv'))
field <- read_csv(paste0(data, 'field_experiment_cleaned.csv'))
field_raw <- read_csv(paste0(data, '/field_raw_trials.csv'))
First, how many participants?
nrow(web)
## [1] 903
Sex division
table(web$Sex)
##
## F M
## 681 222
Ages
summary(web$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 23.00 29.00 32.92 40.00 84.00
Order division
# Counts
table(web$Order)
##
## l_first r_first
## 436 467
# Percentage
prop.table(table(web$Order)) * 100
##
## l_first r_first
## 48.2835 51.7165
First, how many languages?
web %>% count(Name) %>% nrow()
## [1] 25
Does this number correspond with the L1s?
web %>% count(Language) %>% nrow()
## [1] 25
How many families?
web %>% count(Family) %>% nrow()
## [1] 9
How many have the R/L distinction in the L1 among the languages?
web %>% count(Language, r_l_distinction_L1) %>% count(r_l_distinction_L1)
How many really use the alveolar trill in L1 among the languages?
web %>% count(Language, trill_real_L1) %>% count(trill_real_L1)
How many really have the alveolar trill in L1 as an allophone among the languages?
web %>% count(Language, trill_occ_L1) %>% count(trill_occ_L1)
What about the same questions for L2. But this will not neatly sum up to 25, due to various possible scenarios for L2 within a specific L1.
How many have the R/L distinction in the L2 among the languages?
web %>% count(Language, r_l_distinction_L2) %>% count(r_l_distinction_L2)
How many really use the alveolar trill in L2 among the languages?
web %>% count(Language, trill_real_L2) %>% count(trill_real_L2)
How many really have the alveolar trill in L2 as an allophone among the languages?
web %>% count(Language, trill_occ_L2) %>% count(trill_occ_L2)
What is the grand average congruent behavior?
mean(web$Match)
## [1] 0.8726467
87.3%
What about only among those who have L1 without the distinction?
web %>%
filter(r_l_distinction_L1 == "0") %>%
summarize(mean_match = mean(Match, na.rm = TRUE))
83.9%
What about only among those who have L1 without the distinction and no L2 that distinguishes?
web %>%
filter(r_l_distinction_L1 == "0") %>%
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L1 == '0') %>%
summarize(mean_match = mean(Match, na.rm = TRUE))
85.1%
Compute average matching behavior per language and sort:
web_avg <- web %>%
group_by(Language) %>%
summarize(M = mean(Match)) %>%
arrange(desc(M)) %>%
mutate(percent = round(M, 2) * 100,
percent = str_c(percent, '%'))
# Show:
web_avg %>% print(n = Inf)
## # A tibble: 25 × 3
## Language M percent
## <chr> <dbl> <chr>
## 1 FI 1 100%
## 2 FR 0.982 98%
## 3 EN 0.974 97%
## 4 DE 0.953 95%
## 5 SE 0.952 95%
## 6 DK 0.944 94%
## 7 HU 0.941 94%
## 8 JP 0.927 93%
## 9 KR 0.909 91%
## 10 GR 0.905 90%
## 11 PL 0.887 89%
## 12 RU 0.872 87%
## 13 EE 0.860 86%
## 14 FA 0.857 86%
## 15 AM 0.85 85%
## 16 ZU 0.85 85%
## 17 IT 0.846 85%
## 18 TR 0.811 81%
## 19 ES 0.806 81%
## 20 GE 0.8 80%
## 21 TH 0.8 80%
## 22 PT 0.770 77%
## 23 RO 0.742 74%
## 24 AL 0.7 70%
## 25 CN 0.696 70%
Check some demographics, also to report in the paper. First, the number of participants per language:
web %>%
count(Name, sort = TRUE) %>% print(n = Inf)
## # A tibble: 25 × 2
## Name n
## <chr> <int>
## 1 German 85
## 2 Portuguese 61
## 3 French 57
## 4 Japanese 55
## 5 Polish 53
## 6 Italian 52
## 7 Russian 47
## 8 Chinese 46
## 9 Estonian 43
## 10 Greek 42
## 11 English 39
## 12 Turkish 37
## 13 Spanish 36
## 14 Hungarian 34
## 15 Romanian 31
## 16 Korean 22
## 17 Farsi 21
## 18 Swedish 21
## 19 Armenian 20
## 20 Thai 20
## 21 Zulu 20
## 22 Danish 18
## 23 Finnish 18
## 24 Georgian 15
## 25 Albanian 10
Then, the number of L1 speakers who have R/L distinction vs. who don’t:
web %>% count(r_l_distinction_L1) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
How many people do not have any L2?
sum(is.na(web$L2)) / nrow(web)
## [1] 0.1351052
# bilinguals:
1 - sum(is.na(web$L2)) / nrow(web)
## [1] 0.8648948
Check how many people knew English as their L2:
web %>% count(EnglishL2YesNo) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
Of those that don’t use a R/L distinction in their L1, how many use R/L distinction in their L2? (double-check if logic alright!)
web %>%
filter(r_l_distinction_L1 == 0) %>%
count(r_l_distinction_L2 == 1) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
How many “pure” speakers were there?, i.e., those people that 1) don’t know English, 2) don’t use an L1 with a R/L distinction, and 3) don’t know an L2 that distinguishes R/L.
web %>%
filter(r_l_distinction_L1 == 0) %>% # excludes English as well
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L2 == 0) %>%
nrow()
## [1] 1
1 person!
Let’s check if this is correct. This gives the list of all participants for whom this applies.
web %>%
filter(r_l_distinction_L1 == 0 & !EnglishL2YesNo & r_l_distinction_L2 == 0) %>%
print(n = 50)
## # A tibble: 1 × 18
## ID Match Language Sex Age Name Script Family Autotyp_Area L2
## <chr> <dbl> <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr>
## 1 JP_2040343 1 JP F 48 Japane… diffe… Japan… N Coast Asia chin…
## # ℹ 8 more variables: EnglishL2YesNo <lgl>, Order <chr>,
## # r_l_distinction_L1 <dbl>, trill_real_L1 <dbl>, trill_occ_L1 <dbl>,
## # r_l_distinction_L2 <dbl>, trill_real_L2 <dbl>, trill_occ_L2 <dbl>
Are these really “pure”? What languages do they speak?
web %>%
filter(r_l_distinction_L1 == 0) %>% # excludes English as well
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L2 == 0) %>%
count(Language)
One Japanese speaker.
Nevertheless, let’s explore whether these also show matches?
web %>%
filter(r_l_distinction_L1 == 0) %>% # excludes English as well
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L2 == 0) %>%
count(Match) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
Yes, similar to above 85%.
Now, some relevant plots:
Distribution of matches in the web data by Language.
Distribution of matches in the web data by Family.
Distribution of matches in the web data by Area.
Distribution of matches in the web data by [r]/[l] distinction in L1.
Distribution of matches in the web data by [r] is main allophone in L1.
Distribution of matches in the web data by [r]/[l] distinction in L2.
Distribution of matches in the web data by [r] is main allophone in L2.
Distribution of matches in the web data (i.e., not modelled but the actual raw data) by Order (columns) and the [r] is main allophone in the L1 (rows).
Our experiment has two trials and counterbalancing of order, but we code only strict matches, meaning that there are four combination of responses for two values of Match: 1/1 = match, 1/0, 0/1 = partial matches (i.e., no match), and 0/0 = mismatch (i.e., no match). If we assume that that the choices for the two sounds are independent, then We can simulate this:
sim_chance_level <- replicate(1000, # simulate 1000 "experiments"
mean(replicate(1000, # the proportion of "matches" for 1000 random "participants"
("r" == sample(c("r","l"), size=1, prob=c(0.5, 0.5))) &&
("l" == sample(c("r","l"), size=1, prob=c(0.5, 0.5)))))); # only full matches matter (here order does not matter as the repose is random
hist(100*sim_chance_level, main="Histogram of chance level", xlab="% full matches", col="gray90", xlim=c(0,100));
abline(v=100*mean(sim_chance_level), col="blue", lwd=2);
Figure 1. Simulating the % of ‘full’
matches across 1000 ‘experiments’ each with 1000 ‘participants’
that associate randomly (at 50%) [r] and [l] with the wavy or the
straight line, assuming independence between the two choices. Figure
generated using R
version 4.3.3 (2024-02-29)
Indeed, as expected, the chance level of such a ‘full’, conservative match is indeed 25%.
However, it is clear that we cannot assume independence, as the choice made for the first sound clearly affects the choice made for the second, but it is unclear how much.
At the other extreme, if we assume perfect dependence between the two choices (i.e., whatever is chosen for the first sound, will not be chosen for the second), we have the following baseline probability:
sim_chance_level2 <- replicate(1000, # simulate 1000 "experiments"
mean(replicate(1000, # the proportion of "matches" for 1000 random "participants"
("r" == sample(c("r","l"), size=1, prob=c(0.5, 0.5)))))); # only the first choice really matters (the second is fully dependent on it, so if the first is correct, the second will be as well, resulting in a match, and vice-versa)
hist(100*sim_chance_level2, main="Histogram of chance level", xlab="% full matches", col="gray90", xlim=c(0,100));
abline(v=100*mean(sim_chance_level2), col="blue", lwd=2);
Figure 2. Simulating the % of ‘full’
matches across 1000 ‘experiments’ each with 1000 ‘participants’
that associate randomly (at 50%) [r] and [l] with the wavy or the
straight line, assuming perfect dependence between the two choices.
Figure generated using R version 4.3.3
(2024-02-29)
Indeed, as expected, the chance level of such a ‘full’, conservative match is indeed 50%.
It is clear that we are in between these extremes, probably closer to 50% than to 25%, so we will use the 50% chance level throughout being very clear that this is a conservative level.
([)Please note that we abstain from attempting to estimating the chance level from the data, as it is very hard to disentangle, in our experimental design, the chance level from the sound symbolism bias we want to estimate.)
In the various logistic regression models that we fit here, we use various techniques for judging the contribution of any given fixed effect:
bayestestR::ci with method = "ETI" for
details), which should not, in principle, include 0,brms::hypothesis
for details): these are always directional and in the same direction as
the point estimate, and should not be taken, except in the case of the
intercept, as being a priori justified, but simply as a way of
summarizing where the bulk of the posterior distribution is relative to
0; we use the default cut-off of 0.95 to judge the “existence” of the
effect in the given direction (marked with a star “*“).Here we are interested in predicting the probability of a perfect match, Match, from a set of potential predictors and their interactions. A priori, given our hypothesis, the potential predictors here are the order of presentation, Order, and the various properties of the L1(s) and L2(s) spoken by the participant in what concerns [r] and [l], namely:
Importantly, we have no a priori reason to consider the participant’s Sex or Age as potential predictors, but we can nevertheless test their influence on Match.
We code here all the [r] and [r]/[l] predictors as factors with a treatment contrast as follows:
On the other hand, we code Order as a factor with possible values “r_first” and “l_first” but using a contrast (or deviation) coding (-0.5, +0.5), because Order is (almost) balanced at 51.7% “r_first” vs 48.3% “l_first”, which ensures that the intercept represents (roughly) the grand mean (see, for example, Chapter Sum contrasts in An Introduction to Bayesian Data Analysis for Cognitive Science by Bruno Nicenboim, Daniel Schad, and Shravan Vasishth). (Please note that the other predictors are not balanced, making the treatment contrast better suited.)
# Make sure we code the factors with the intended baseline levels and contrasts:
web$Order <- factor(web$Order, levels = c("r_first", "l_first")); contrasts(web$Order) <- c(-0.5, +0.5)
web$r_l_distinction_L1 <- factor(c("same", "distinct")[web$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
web$trill_real_L1 <- factor(c("no", "yes")[web$trill_real_L1 + 1], levels=c("no", "yes"));
web$trill_occ_L1 <- factor(c("no", "yes")[web$trill_occ_L1 + 1], levels=c("no", "yes"));
web$r_l_distinction_L2 <- factor(c("same", "distinct")[web$r_l_distinction_L2 + 1], levels=c("same", "distinct"));
web$trill_real_L2 <- factor(c("no", "yes")[web$trill_real_L2 + 1], levels=c("no", "yes"));
web$trill_occ_L2 <- factor(c("no", "yes")[web$trill_occ_L2 + 1], levels=c("no", "yes"));
web %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# match is quite unbalanced: 12.7% vs 87.3%
web %>% count(Order) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# the order effect is decently balanced: 51.6% vs 48.3%
web %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# highly imbalanced: 15.8% vs 84.2%
web %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# ok: 58.8% vs 41.2%
web %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 100% are 1 --> excluded from the model
## And for L2, just in case
web %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 13.5% missing, the rest almost all (86.3%) are 1 and 0.2% are 0 ---> exclude as well
web %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 13.5% missing, the rest is balanced: 53.8% vs 32.7%
web %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 13.5% missing, the rest are all (86.5%) are 1 ---> exclude as well
It can be seen that, in our data, trill_occ_L1 and trill_occ_L2 are both 100% “yes” (for the non-missing data cases), which means that we will not include them in any further analyses.
Because we want to predict the probability of a Match, which is a binary variable, we will use logistic regression throughout. However, logistic regression coefficients are notoriously hard to interpret, as they are reported on the log odds scale. Here, we show the intercepts and regression slopes both as log odds ratios and as probabilities, as appropriate. (See, for example, here, here and here for more explanations.)
Please note that for the intercept (α) this represents the probability of a match when all the predictors are at 0.0 or their baseline levels; for example, an intercept of 2.87 corresponds to a probability of 94.6% of a match when holding all the other predictors at 0 or at base level.
For the slopes (β), this interpretation changes, and we show both the equivalent percent change in odds and the change in the probability of a match between this predictor’s baseline and non-baseline levels (or for a unit increase from 0.0 to 1.0) when all the other predictors are at baseline (or at 0.0). For example, a slope of -0.84 for a binary DV corresponds to an odds ratio of 0.43, which represents a decrease of 56.8% in the odds of a match, or, equivalently, a decrease by 6.2% in the probability of a match from the baseline 94.6% when the DV is not at its baseline level to 88.4%.
An important point to clarify is what should be the random effects structure of our regression models.
A priori, there are three potential “grouping factors”: Language, Family and (AUTOTYP) Area, that are meaningfully considered as random effects. We have 25 L1 languages, from 9 families and 6 areas:
unique(web[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)
However, when it comes to the r/l distinction and the [r] is main allophone in the language (L1 or L2), it is clear that these variables do not, by definition, vary within Language (so there should be no random slope here), and, at least in our current data, they very very little within the levels of the other two grouping factors as well (see below), raising the legitimate question whether we should model random slopes for these two variables at all.
Order: we know this varies within all levels of Language, Family and Area because of the experimental design, so Order should have varying slopes for all three.
r/l distinction in L1 (r_l_distinction_L1):
table(web$r_l_distinction_L1, web$Language)
##
## AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH
## same 0 0 46 0 0 0 0 0 0 0 0 0 0 0 0 55 22 0 0 0 0 0 0
## distinct 10 20 0 85 18 43 39 36 21 18 57 15 42 34 52 0 0 53 61 31 47 21 20
##
## TR ZU
## same 0 20
## distinct 37 0
table(web$r_l_distinction_L1, web$Family)
##
## Benue-Congo Finno-Ugric IE Japanese Kartvelian Korean Sino-Tibetan
## same 20 0 0 55 0 22 46
## distinct 0 95 593 0 15 0 0
##
## Tai-Kadai Turkic
## same 0 0
## distinct 20 37
table(web$r_l_distinction_L1, web$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
## same 0 0 0 77 20
## distinct 539 93 108 0 0
##
## Southeast Asia
## same 46
## distinct 20
This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.
r/l distinction in L2 (r_l_distinction_L2):
table(web$r_l_distinction_L2, web$Language)
##
## AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH
## same 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
## distinct 9 20 30 84 16 43 28 34 18 18 46 15 42 26 50 28 21 52 42 30 43 19 18
##
## TR ZU
## same 0 0
## distinct 29 18
table(web$r_l_distinction_L2, web$Family)
##
## Benue-Congo Finno-Ugric IE Japanese Kartvelian Korean Sino-Tibetan
## same 0 0 1 1 0 0 0
## distinct 18 87 533 28 15 21 30
##
## Tai-Kadai Turkic
## same 0 0
## distinct 18 29
table(web$r_l_distinction_L2, web$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
## same 1 0 0 1 0
## distinct 478 82 104 49 18
##
## Southeast Asia
## same 0
## distinct 48
There is almost no “same” at all, so random slopes are arguably not justified.
[r] is main allophone in L1 (trill_real_L1):
table(web$trill_real_L1, web$Language)
##
## AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH TR
## no 0 0 46 82 18 0 38 0 0 0 56 14 42 0 0 55 22 0 61 0 0 20 20 37
## yes 10 20 0 3 0 43 1 36 21 18 1 1 0 34 52 0 0 53 0 31 47 1 0 0
##
## ZU
## no 20
## yes 0
table(web$trill_real_L1, web$Family)
##
## Benue-Congo Finno-Ugric IE Japanese Kartvelian Korean Sino-Tibetan
## no 20 0 317 55 14 22 46
## yes 0 95 276 0 1 0 0
##
## Tai-Kadai Turkic
## no 20 37
## yes 0 0
table(web$trill_real_L1, web$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
## no 317 51 0 77 20
## yes 222 42 108 0 0
##
## Southeast Asia
## no 66
## yes 0
This is almost uniform within Language, but there is variation for the IE level of Family (almost 50%:50% between “no” and “yes”), and between 2 levels of Area (Europe with 60%:40% and Greater Mesopotamia with 55%:45% “no”:“yes”), but this is not enough to justify random slopes either.
[r] is main allophone in L2 (trill_real_L2):
table(web$trill_real_L2, web$Language)
##
## AL AM CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PL PT RO RU SE TH TR
## no 5 1 29 38 13 5 16 22 16 8 28 3 28 18 30 26 21 34 19 21 32 11 18 28
## yes 4 19 1 46 3 38 13 12 2 10 18 12 14 8 20 3 0 18 23 9 11 8 0 1
##
## ZU
## no 16
## yes 2
table(web$trill_real_L2, web$Family)
##
## Benue-Congo Finno-Ugric IE Japanese Kartvelian Korean Sino-Tibetan
## no 16 31 314 26 3 21 29
## yes 2 56 220 3 12 0 1
##
## Tai-Kadai Turkic
## no 18 28
## yes 0 1
table(web$trill_real_L2, web$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Coast Asia S Africa
## no 283 48 45 47 16
## yes 196 34 59 3 2
##
## Southeast Asia
## no 47
## yes 1
Here there is enough variation within the levels for all three factors, justifying the inclusion of random slopes.
Given these, we did the following:
if( !file.exists("./cached_results/b1_res_web.RData") )
{
# prior predictive checks:
# what priors we can set (use Order for that):
get_prior(Match ~ 1 + Order +
(1 + Order | Language) +
(1 + Order | Family) +
(1 + Order | Autotyp_Area),
data=web,
family=bernoulli(link='logit'));
# -> "sd" priors seem alright (student_t(3, 0, 2.5)), as does the "Intercept" (student_t(3, 0, 2.5)), but lkj(1) for "cor" might too accepting of extreme correlations, and (flat) for "b" is clearly not ok
# so, we keep "sd" and "Intercept" but use lkj(2) for "cor" and student_t(5, 0, 2.5) for "b"
b_priors <- brms::brm(Match ~ 1 + Order +
(1 + Order | Language) +
(1 + Order | Family) +
(1 + Order | Autotyp_Area),
data=web,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"),
brms::set_prior("lkj(2)", class="cor")),
sample_prior='only', # needed for prior predictive checks
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is a bit extreme
# the model:
b1 <- brm(Match ~ 1 + Order +
(1 + Order | Language) +
(1 + Order | Family) +
(1 + Order | Autotyp_Area),
data=web,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"),
brms::set_prior("lkj(2)", class="cor")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
mcmc_plot(b1, type="trace");
bayestestR::hdi(b1, ci=0.95);
hypothesis(b1, "Order1 = 0");
# posterior predictive checks
pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
ncol=1); # seems fine (the observed, y, does fall in the predicted distributions, y_{rep} for the max, mean and min)
# save the model:
b1_res_web <- list("b_priors"=b_priors, "b1"=b1);
saveRDS(b1_res_web, "./cached_results/b1_res_web.RData", compress="xz");
} else
{
b1_res_web <- readRDS("./cached_results/b1_res_web.RData");
}
Priors As per various guidelines (see, for example,
https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations)
and experiments, we will use student_t(5, 0, 2.5) for the
slopes β, and lkj(2) for the random effects
correlation structure, keeping the default brms priors for
everything else, as they seem both sane and to perform well in our prior
predictive checks. So, the priors we use are:
get_prior(b1_res_web$b_priors);
grid.arrange(pp_check(b1_res_web$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b1_res_web$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b1_res_web$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is a bit extreme
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 3. Prior predictive checks. Figure generated
using R version
4.3.3 (2024-02-29)
and it can be seen that they are quite ok, even if the observed p(match) is a bit extreme, but still within what the prior distribution covers.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b1_res_web$b1);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + Order + (1 + Order | Language) + (1 + Order | Family) + (1 + Order | Autotyp_Area)
## Data: web (Number of observations: 903)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.43 0.38 0.02 1.40 1.00 7044
## sd(Order1) 0.68 0.62 0.02 2.31 1.00 8744
## cor(Intercept,Order1) 0.05 0.45 -0.80 0.84 1.00 10857
## Tail_ESS
## sd(Intercept) 9633
## sd(Order1) 10223
## cor(Intercept,Order1) 10571
##
## ~Family (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.36 0.30 0.01 1.12 1.00 8623
## sd(Order1) 0.69 0.55 0.03 2.06 1.00 8960
## cor(Intercept,Order1) 0.03 0.45 -0.80 0.83 1.00 10376
## Tail_ESS
## sd(Intercept) 9633
## sd(Order1) 9252
## cor(Intercept,Order1) 10974
##
## ~Language (Number of levels: 25)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.55 0.21 0.16 1.02 1.00 6419
## sd(Order1) 0.93 0.36 0.26 1.70 1.00 7836
## cor(Intercept,Order1) 0.42 0.33 -0.31 0.92 1.00 8766
## Tail_ESS
## sd(Intercept) 5539
## sd(Order1) 6938
## cor(Intercept,Order1) 9511
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.03 0.36 1.28 2.73 1.00 9372 9214
## Order1 -1.04 0.61 -2.34 0.12 1.00 10071 10084
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b1_res_web$b1, type="trace");
## No divergences to plot.
Figure 4. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks seems to be perfectly fine:
pp_check(b1_res_web$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 5. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(pp_check(b1_res_web$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
pp_check(b1_res_web$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
pp_check(b1_res_web$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 6. Posterior predictive distribution. Figure
generated using R
version 4.3.3 (2024-02-29)
Fixed effects:
mcmc_plot(b1_res_web$b1);
Figure 7. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
plot_model(b1_res_web$b1, type="emm", terms=c("Order"));
Figure 8. Fixed effects. Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 2.03 95%CI [1.28, 2.73] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=88.4% with 95%CI [78.3%, 93.9%] ≫ 50%.
Order The slope of Order is -1.04 95%CI [-2.34, 0.12] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.96*. In terms of probabilities, presenting [l] first (Order == “l_first”) results in a decrease by 15.5% [0.666%, 52.4%] in the probability of a match from the 88.4% [78.3%, 93.9%] when [r] is presented first 73.0% [38.3%, 91.9%] when presenting [r] first.
Random effects:
#get_variables(b1_res_web$b1);
grid.arrange(
# 1 | Language:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Family:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Autotyp_Area:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# Order | Language:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
theme(legend.position="bottom") +
NULL,
# Order | Family:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
# Order | Autotyp_Area:
b1_res_web$b1 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
nrow=2);
Figure 9. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match in the online experiment, while controlling for Language, Family and Area, is clearly much higher than the conservative change level of 50%, being estimated as about 88% with a 95%CI of about [78%, 94%], showing that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line. The order of presentation does seem to matter, in the sense that presenting [r] first increases the chance of a match by about 15% with a wide 95% [0.7%, 52%].
if( !file.exists("./cached_results/b2_res_web.RData") )
{
# the model:
b2 <- brm(Match ~ 1 + Order + # random slope
r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + # only random intercepts for these
trill_real_L2 + # but random slope for this one
(1 + Order + trill_real_L2 | Language) +
(1 + Order + trill_real_L2 | Family) +
(1 + Order + trill_real_L2 | Autotyp_Area),
data=web,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"),
brms::set_prior("lkj(2)", class="cor")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
mcmc_plot(b2, type="trace");
bayestestR::hdi(b2, ci=0.95);
hypothesis(b2, c("Order1 = 0", "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0"));
# posterior predictive checks
pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1); # seems perfect
# save the model:
b2_res_web <- list("b2"=b2);
saveRDS(b2_res_web, "./cached_results/b2_res_web.RData", compress="xz");
} else
{
b2_res_web <- readRDS("./cached_results/b2_res_web.RData");
}
Priors are as before.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b2_res_web$b2);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + Order + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 + Order + trill_real_L2 | Language) + (1 + Order + trill_real_L2 | Family) + (1 + Order + trill_real_L2 | Autotyp_Area)
## Data: web (Number of observations: 781)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.60 0.55 0.02 2.00 1.00
## sd(Order1) 1.02 0.88 0.04 3.28 1.00
## sd(trill_real_L2yes) 1.16 1.00 0.05 3.74 1.00
## cor(Intercept,Order1) 0.08 0.42 -0.72 0.81 1.00
## cor(Intercept,trill_real_L2yes) -0.02 0.41 -0.77 0.74 1.00
## cor(Order1,trill_real_L2yes) -0.00 0.42 -0.76 0.76 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 7159 9073
## sd(Order1) 6980 8396
## sd(trill_real_L2yes) 7534 8426
## cor(Intercept,Order1) 10077 10771
## cor(Intercept,trill_real_L2yes) 10390 10936
## cor(Order1,trill_real_L2yes) 10479 10735
##
## ~Family (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.86 0.58 0.04 2.23 1.00
## sd(Order1) 1.05 0.79 0.04 2.98 1.00
## sd(trill_real_L2yes) 0.91 0.85 0.03 3.13 1.00
## cor(Intercept,Order1) 0.14 0.41 -0.68 0.83 1.00
## cor(Intercept,trill_real_L2yes) -0.08 0.41 -0.81 0.71 1.00
## cor(Order1,trill_real_L2yes) -0.03 0.41 -0.78 0.75 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5893 7042
## sd(Order1) 6711 7431
## sd(trill_real_L2yes) 9535 9792
## cor(Intercept,Order1) 9073 9221
## cor(Intercept,trill_real_L2yes) 10272 9746
## cor(Order1,trill_real_L2yes) 10743 10102
##
## ~Language (Number of levels: 25)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.39 0.27 0.02 1.01 1.00
## sd(Order1) 0.91 0.46 0.09 1.88 1.00
## sd(trill_real_L2yes) 0.79 0.47 0.05 1.85 1.00
## cor(Intercept,Order1) 0.21 0.39 -0.63 0.85 1.00
## cor(Intercept,trill_real_L2yes) -0.07 0.40 -0.77 0.71 1.00
## cor(Order1,trill_real_L2yes) 0.19 0.37 -0.58 0.82 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5511 7578
## sd(Order1) 5891 4854
## sd(trill_real_L2yes) 5811 5807
## cor(Intercept,Order1) 5995 8490
## cor(Intercept,trill_real_L2yes) 8223 9753
## cor(Order1,trill_real_L2yes) 8951 10113
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 3.58 2.63 -0.67 9.41 1.00 9134
## Order1 -0.87 0.80 -2.54 0.68 1.00 9786
## r_l_distinction_L1distinct 0.25 0.95 -1.71 2.05 1.00 9300
## r_l_distinction_L2distinct -1.16 2.53 -6.80 2.83 1.00 8942
## trill_real_L1yes -1.05 0.42 -1.91 -0.25 1.00 10194
## trill_real_L2yes 0.13 0.95 -1.61 2.26 1.00 8981
## Tail_ESS
## Intercept 7826
## Order1 9854
## r_l_distinction_L1distinct 10079
## r_l_distinction_L2distinct 7859
## trill_real_L1yes 10650
## trill_real_L2yes 8416
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b2_res_web$b2, type="trace");
## No divergences to plot.
Figure 10. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks seems to be perfectly fine:
pp_check(b2_res_web$b2, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 11. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(pp_check(b2_res_web$b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2_res_web$b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2_res_web$b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 12. Posterior predictive distribution (the
x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3
(2024-02-29)
Fixed effects:
mcmc_plot(b2_res_web$b2);
Figure 13. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(plot_model(b2_res_web$b2, type="emm", terms=c("Order")),
plot_model(b2_res_web$b2, type="emm", terms=c("r_l_distinction_L1")),
plot_model(b2_res_web$b2, type="emm", terms=c("r_l_distinction_L2")),
plot_model(b2_res_web$b2, type="emm", terms=c("trill_real_L1")),
plot_model(b2_res_web$b2, type="emm", terms=c("trill_real_L2")),
ncol=2);
Figure 14. Fixed effects. Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 3.58 95%CI [-0.67, 9.41] on the log odds ratio (LOR) scale, which is > 0 but has very large uncertainties so that the 95%CI includes 0; formal hypothesis testing against 0 is p(α>0)=0.95; this translates into a probability of a match p(match)=97.3% 95%CI [33.9%, 100.0%] which includes 50% in the 95%CI.
Order The slope of Order is -0.87 95%CI [-2.54, 0.68] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.89. In terms of probabilities, presenting [l] first (Order1 == “l_first”) results in a decrease by 3.5% [0.004%, 30.0%] in the probability of a match from the 97.3% [33.9%, 100.0%] when [r] is presented first 93.8% [14.8%, 100.0%] when [l] is presented first.
Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is 0.25 95%CI [-1.71, 2.05] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.62.
Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -1.16 95%CI [-6.80, 2.83] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.66.
Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -1.05 95%CI [-1.91, -0.25] on the log odds ratio (LOR) scale, which is clearly < 0; formal hypothesis testing against 0 is p(β<0)=0.99*. In terms of probabilities, having [r] as the main allophone (trill_real_L1 == “yes”) results in a decrease by 4.6% [0.002%, 26.8%] in the probability of a match from the 97.3% [33.9%, 100.0%] when [r] is not the primary allophone to 92.6% [14.6%, 100.0%] when it is.
Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is 0.13 95%CI [-1.61, 2.26] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.53.
Random effects:
#get_variables(b2_res_web$b2);
grid.arrange(
# 1 | Language:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Family:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Autotyp_Area:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# Order | Language:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
theme(legend.position="bottom") +
NULL,
# Order | Family:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
# Order | Autotyp_Area:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Order1") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("Order | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Language:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Language[Language, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Family:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Family[Family, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Autotyp_Area:
b2_res_web$b2 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
nrow=3);
Figure 15. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match in the online experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, Family and Area, seem to still be higher than the conservative change level of 50%, being estimated as about 97% but with a much wider 95%CI of about [34%, 100.0%] which now does include 50%. The only predictor that seems to make a difference is if [r] is the main allophone in the L1(s) known by the participant, in that those participants that have an L1 with [r] as the main allophone have a slightly smaller probability of a perfect match by about 4.5% [0%, 27%] to about 93% [15%, 100%] than the others. Here, order does not seem to make clear a difference, but there is still a suggestion that presenting [r] first increases the probability of a match by about 3.5% [0%, 30%].
First, how many participants?
nrow(field)
## [1] 127
Sex division
table(field$Sex)
##
## f m
## 91 36
Ages
summary(field$Age)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 19.00 20.00 28.55 34.50 75.00
First, how many languages?
field %>% count(Language) %>% nrow()
## [1] 6
Does this number correspond with the L1s?
field %>% count(Name) %>% nrow()
## [1] 6
How many families?
field %>% count(Family) %>% nrow()
## [1] 4
How many have the R/L distinction in the L1 among the languages?
field %>% count(Name, r_l_distinction_L1) %>% count(r_l_distinction_L1)
How many really use the alveolar trill in L1 among the languages?
field %>% count(Name, trill_real_L1) %>% count(trill_real_L1)
How many really have the alveolar trill in L1 as an allophone among the languages?
field %>% count(Name, trill_occ_L1) %>% count(trill_occ_L1)
What about the same questions for L2. But this will not neatly sum up to 25, due to various possible scenarios for L2 within a specific L1.
How many have the R/L distinction in the L2 among the languages?
field %>% count(Name, r_l_distinction_L2) %>% count(r_l_distinction_L2)
How many really use the alveolar trill in L2 among the languages?
field %>% count(Name, trill_real_L2) %>% count(trill_real_L2)
How many really have the alveolar trill in L2 as an allophone among the languages?
field %>% count(Name, trill_occ_L2) %>% count(trill_occ_L2)
What is the grand average congruent behavior?
mean(field$Match)
## [1] 0.976378
97%!!!
What about only among those who have L1 without the distinction?
field %>%
filter(r_l_distinction_L1 == "0") %>%
summarize(mean_match = mean(Match, na.rm = TRUE))
100%! WOW.
What about only among those who have L1 without the distinction and no L2 that distinguishes?
field %>%
filter(r_l_distinction_L1 == "0") %>%
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L2 == '0') %>%
summarize(mean_match = mean(Match, na.rm = TRUE))
There are no such people.
Compute average matching behavior per language and sort:
field_avg <- field %>%
group_by(Language) %>%
summarize(M = mean(Match)) %>%
arrange(desc(M)) %>%
mutate(percent = round(M, 2) * 100,
percent = str_c(percent, '%'))
# Show:
field_avg %>% print(n = Inf)
## # A tibble: 6 × 3
## Language M percent
## <chr> <dbl> <chr>
## 1 BR 1 100%
## 2 PA 1 100%
## 3 VA 1 100%
## 4 BE 0.982 98%
## 5 DE 0.947 95%
## 6 SR 0.923 92%
Check some demographics, also to report in the paper. First, the number of participants per language:
field %>%
count(Name, sort = TRUE) %>% print(n = Inf)
## # A tibble: 6 × 2
## Name n
## <chr> <int>
## 1 English UK 55
## 2 Tashlhiyt Berber 20
## 3 German 19
## 4 Brazilian Portuguese 13
## 5 Daakie 12
## 6 Palikur 8
Then, the number of L1 speakers who have R/L distinction vs. who don’t:
field %>% count(r_l_distinction_L1) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
How many people do not have any L2?
# raw count of no L2; raw count of all ppl
sum(is.na(field$L2)); nrow(field)
## [1] 52
## [1] 127
# percentage no L2
sum(is.na(field$L2)) / nrow(field)
## [1] 0.4094488
# percentage with L2
1 - sum(is.na(field$L2)) / nrow(field)
## [1] 0.5905512
Check how many people knew English as their L2:
field %>% count(EnglishL2YesNo) %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
Of those that don’t use a R/L distinction in their L1, how many use R/L distinction in their L2? (double-check if logic alright!)
field %>%
filter(r_l_distinction_L1 == '0') %>%
count(r_l_distinction_L2 == '1') %>%
mutate(prop = n / sum(n),
percent = round(prop, 2) * 100,
percent = str_c(percent, '%'))
How many “pure” speakers were there?, i.e., those people that 1) don’t know English, 2) don’t use an L1 with a R/L distinction, and 3) don’t know an L2 that distinguishes R/L.
field %>%
filter(r_l_distinction_L1 == '0') %>% # excludes English as well
filter(!EnglishL2YesNo) %>%
filter(r_l_distinction_L2 == '0') %>%
nrow()
## [1] 0
None.
Now, some relevant plots:
Distribution of matches in the field data by Language.
Distribution of matches in the field data by Family.
Distribution of matches in the field data by Area.
Distribution of matches in the field data by [r]/[l] distinction in L1.
Distribution of matches in the field data by [r] is main allophone in L1.
Distribution of matches in the field data by [r]/[l] distinction in L2.
Distribution of matches in the field data by [r] is main allophone in L2.
We use the same ideas as for the web experiment above, except that now Order does not exist anymore.
# Make sure we code the factors with the intended baseline levels and contrasts:
field$r_l_distinction_L1 <- factor(c("same", "distinct")[field$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
field$trill_real_L1 <- factor(c("no", "yes")[field$trill_real_L1 + 1], levels=c("no", "yes"));
field$trill_occ_L1 <- factor(c("no", "yes")[field$trill_occ_L1 + 1], levels=c("no", "yes"));
field$r_l_distinction_L2 <- factor(c("same", "distinct")[field$r_l_distinction_L2 + 1], levels=c("same", "distinct"));
field$trill_real_L2 <- factor(c("no", "yes")[field$trill_real_L2 + 1], levels=c("no", "yes"));
field$trill_occ_L2 <- factor(c("no", "yes")[field$trill_occ_L2 + 1], levels=c("no", "yes"));
field %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# match is quite unbalanced: 2.4% vs 97.6%
field %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# highly imbalanced: 6.3% vs 93.7%
field %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# imbalanced: 74.8% vs 25.2%
field %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# imbalanced: 6.3% vs 93.7%
## And for L2, just in case
field %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 41.7% missing, the rest: 0.8% vs 57.5%
field %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 41.7% missing, the rest: 22% vs 36.2%
field %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
# 41.7% missing, the rest are all (58.3%) are "yes" ---> exclude
It can be seen that, in our data, trill_occ_L2 is 100% “yes” (for the non-missing data cases), which means that we will not include it in any further analyses.
However, r_l_distinction_L1 and trill_occ_L1 are perfectly correlated:
table(field$r_l_distinction_L1, field$trill_occ_L1);
##
## no yes
## same 8 0
## distinct 0 119
so we will only consider r_l_distinction_L1.
As for the web experiment above:
unique(field[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)
It seems that including Family and Area do not make any sense in this case.
r/l distinction in L1 (r_l_distinction_L1):
table(field$r_l_distinction_L1, field$Language)
##
## BE BR DE PA SR VA
## same 0 0 0 8 0 0
## distinct 55 20 19 0 13 12
table(field$r_l_distinction_L1, field$Family)
##
## Afro-Asiatic Arawakan Austronesian IE
## same 0 8 0 0
## distinct 20 0 12 87
table(field$r_l_distinction_L1, field$Autotyp_Area)
##
## Europe N Africa NE South America Oceania
## same 0 0 8 0
## distinct 87 20 0 12
This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.
r/l distinction in L2 (r_l_distinction_L2):
table(field$r_l_distinction_L2, field$Language)
##
## BE BR DE PA SR VA
## same 1 0 0 0 0 0
## distinct 18 20 14 8 1 12
table(field$r_l_distinction_L2, field$Family)
##
## Afro-Asiatic Arawakan Austronesian IE
## same 0 0 0 1
## distinct 20 8 12 33
table(field$r_l_distinction_L2, field$Autotyp_Area)
##
## Europe N Africa NE South America Oceania
## same 1 0 0 0
## distinct 33 20 8 12
There is almost no “same” at all, so random slopes are arguably not justified.
[r] is main allophone in L1 (trill_real_L1):
table(field$trill_real_L1, field$Language)
##
## BE BR DE PA SR VA
## no 55 0 19 8 13 0
## yes 0 20 0 0 0 12
table(field$trill_real_L1, field$Family)
##
## Afro-Asiatic Arawakan Austronesian IE
## no 0 8 0 87
## yes 20 0 12 0
table(field$trill_real_L1, field$Autotyp_Area)
##
## Europe N Africa NE South America Oceania
## no 87 0 8 0
## yes 0 20 0 12
This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.
[r] is main allophone in L2 (trill_real_L2):
table(field$trill_real_L2, field$Language)
##
## BE BR DE PA SR VA
## no 9 0 10 8 1 0
## yes 10 20 4 0 0 12
table(field$trill_real_L2, field$Family)
##
## Afro-Asiatic Arawakan Austronesian IE
## no 0 8 0 20
## yes 20 0 12 14
table(field$trill_real_L2, field$Autotyp_Area)
##
## Europe N Africa NE South America Oceania
## no 20 0 8 0
## yes 14 20 0 12
There is variation only within Indo-European/Europe, so random slopes are arguably not justified.
[r] is present in L1 (trill_occ_L1):
table(field$trill_occ_L1, field$Language)
##
## BE BR DE PA SR VA
## no 0 0 0 8 0 0
## yes 55 20 19 0 13 12
table(field$trill_occ_L1, field$Family)
##
## Afro-Asiatic Arawakan Austronesian IE
## no 0 8 0 0
## yes 20 0 12 87
table(field$trill_occ_L1, field$Autotyp_Area)
##
## Europe N Africa NE South America Oceania
## no 0 0 8 0
## yes 87 20 0 12
This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.
Given these, we only include Language as a random effect, and we do not model any random slopes here.
As for the web experiment, but this time without Order (so, the intercept-only model):
if( !file.exists("./cached_results/b1_res_field.RData") )
{
# prior predictive checks:
# what priors we can set:
get_prior(Match ~ 1 +
(1 | Language),
data=field,
family=bernoulli(link='logit'));
b_priors <- brms::brm(Match ~ 1 +
(1 | Language),
data=field,
family=bernoulli(link='logit'),
sample_prior='only', # needed for prior predictive checks
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is quite extreme
# the model:
b1 <- brm(Match ~ 1 +
(1 | Language),
data=field,
family=bernoulli(link='logit'),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
mcmc_plot(b1, type="trace");
bayestestR::hdi(b1, ci=0.95);
# posterior predictive checks
pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
ncol=1); # far from ideal but seems ok...
# save the model:
b1_res_field <- list("b_priors"=b_priors, "b1"=b1);
saveRDS(b1_res_field, "./cached_results/b1_res_field.RData", compress="xz");
} else
{
b1_res_field <- readRDS("./cached_results/b1_res_field.RData");
}
grid.arrange(pp_check(b1_res_field$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b1_res_field$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b1_res_field$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is a bit extreme
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 16. Prior predictive checks. Figure generated
using R version
4.3.3 (2024-02-29)
and it can be seen that they are quite ok, even if the observed p(match) is quite extreme, but still within what the prior distribution covers.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b1_res_field$b1);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + (1 | Language)
## Data: field (Number of observations: 127)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Language (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.98 0.92 0.03 3.42 1.00 7041 7465
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 3.84 0.84 2.37 5.72 1.00 7664 7214
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b1_res_field$b1, type="trace");
## No divergences to plot.
Figure 17. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks seems to be acceptable:
pp_check(b1_res_field$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 18. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(pp_check(b1_res_field$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b1_res_field$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b1_res_field$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 19. Posterior predictive distribution. Figure
generated using R
version 4.3.3 (2024-02-29)
Fixed effects:
mcmc_plot(b1_res_field$b1);
Figure 20. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 3.84 95%CI [2.37, 5.72] on the log odds ratio (LOR) scale, which is clearly » 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=97.9% 95%CI [91.4%, 99.7%] » 50%.
Random effects:
#get_variables(b1_res_field$b1);
grid.arrange(
# 1 | Language:
b1_res_field$b1 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
nrow=1);
Figure 21. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match in the field experiment, while controlling for Language, is much higher than the conservative change level of 50%, being estimated as about 98%, with a 95% credible interval of about [91%, 100%] that excludes 50%, suggesting that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line.
if( !file.exists("./cached_results/b2_res_field.RData") )
{
# the model:
b2 <- brm(Match ~ 1 +
r_l_distinction_L1 + r_l_distinction_L2 +
trill_real_L1 + trill_real_L2 +
(1 | Language),
data=field,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")), # same slope prior as for web
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
mcmc_plot(b2, type="trace");
bayestestR::hdi(b2, ci=0.95);
hypothesis(b2, c("r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0"));
# posterior predictive checks
pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
pp_check(b2, ndraws=100, prefix="ppd") + xlab('p(match)') + ggtitle('Posterior predictive density overlay without the observed data');
grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1); # seems perfect
# save the model:
b2_res_field <- list("b2"=b2);
saveRDS(b2_res_field, "./cached_results/b2_res_field.RData", compress="xz");
} else
{
b2_res_field <- readRDS("./cached_results/b2_res_field.RData");
}
Priors are as before.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b2_res_field$b2);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 | Language)
## Data: field (Number of observations: 74)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Language (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 2.05 2.38 0.06 7.30 1.00 9933 8547
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 10.26 7.57 0.17 27.28 1.00 7762
## r_l_distinction_L1distinct -0.21 2.87 -6.14 5.30 1.00 10425
## r_l_distinction_L2distinct -0.13 3.01 -6.39 5.68 1.00 11392
## trill_real_L1yes -0.13 2.67 -5.49 5.16 1.00 11199
## trill_real_L2yes -0.18 2.65 -5.55 4.98 1.00 10980
## Tail_ESS
## Intercept 5939
## r_l_distinction_L1distinct 9293
## r_l_distinction_L2distinct 9593
## trill_real_L1yes 9418
## trill_real_L2yes 9461
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b2_res_field$b2, type="trace");
## No divergences to plot.
Figure 22. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks are actually pretty good (please note that there are some curves hidden by the observed data, visible in the second plot):
pp_check(b2_res_field$b2, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 23. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
pp_check(b2_res_field$b2, ndraws=500, prefix="ppd") + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 24. Posterior predictive density overlay
without showinf the observed data for n=500 draws. Figure
generated using R
version 4.3.3 (2024-02-29)
grid.arrange(pp_check(b2_res_field$b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2_res_field$b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2_res_field$b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 25. Posterior predictive distribution (the
x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3
(2024-02-29)
Fixed effects:
mcmc_plot(b2_res_field$b2);
Figure 26. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(plot_model(b2_res_field$b2, type="emm", terms=c("r_l_distinction_L1")),
plot_model(b2_res_field$b2, type="emm", terms=c("r_l_distinction_L2")),
plot_model(b2_res_field$b2, type="emm", terms=c("trill_real_L1")),
plot_model(b2_res_field$b2, type="emm", terms=c("trill_real_L2")),
ncol=2);
Figure 27. Fixed effects. Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 10.26 95%CI [0.17, 27.28] on the log odds ratio (LOR) scale, which is > 0 and with a 95%CI that does not includes 0; formal hypothesis testing against 0 is p(α>0)=0.98*; this translates into a probability of a match p(match)=100.0% 95%CI [54.1%, 100.0%] which does not include 50%.
Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is -0.21 95%CI [-6.14, 5.30] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.52.
Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -0.13 95%CI [-6.39, 5.68] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.51.
Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -0.13 95%CI [-5.49, 5.16] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.52.
Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is -0.18 95%CI [-5.55, 4.98] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.53.
Random effects:
#get_variables(b2_res_field$b2);
grid.arrange(
# 1 | Language:
b2_res_field$b2 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
nrow=1);
Figure 28. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match in the field experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, is higher than the conservative chance level of 50%, being estimated as about 100% but with a much wider 95%CI of about [54%, 100.0%] which now includes 50%. No predictors seem to make any difference.
Given that Order does not seem to formally make a difference in the online experiment, we can, in principle, combine the online and the field experiments in a joint analysis by ignoring Order as a predictor (but keeping all the observations).
webfield <- rbind(web[,names(field)], field);
webfield$dataset <- factor(c(rep("web", nrow(web)), rep("field", nrow(field))), levels=c("web", "field"));
webfield$Sex <- toupper(webfield$Sex); # make sex uniform
We use the same ideas as for the web experiment above, except that now Order does not exist anymore.
webfield %>% count(Match) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(r_l_distinction_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(trill_real_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(trill_occ_L1) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(r_l_distinction_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(trill_real_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
webfield %>% count(trill_occ_L2) %>% mutate(prop = paste0(round(100* n / sum(n),1),"%"))
It can be seen that, in our data, trill_occ_L1 and trill_occ_L2 are (almost) at 100% “yes” (for the non-missing data cases), which means that we will not include them in any further analyses.
unique(webfield[,c("Name", "Family", "Autotyp_Area")]) %>% arrange(Autotyp_Area, Family, Name)
r/l distinction in L1 (r_l_distinction_L1):
table(webfield$r_l_distinction_L1, webfield$Language)
##
## AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT
## same 0 0 0 0 46 0 0 0 0 0 0 0 0 0 0 0 0
## distinct 10 20 55 20 0 104 18 43 39 36 21 18 57 15 42 34 52
##
## JP KR PA PL PT RO RU SE SR TH TR VA ZU
## same 55 22 8 0 0 0 0 0 0 0 0 0 20
## distinct 0 0 0 53 61 31 47 21 13 20 37 12 0
table(webfield$r_l_distinction_L1, webfield$Family)
##
## Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric IE
## same 0 8 0 20 0 0
## distinct 20 0 12 0 95 680
##
## Japanese Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
## same 55 0 22 46 0 0
## distinct 0 15 0 0 20 37
table(webfield$r_l_distinction_L1, webfield$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
## same 0 0 0 0 77
## distinct 626 93 108 20 0
##
## NE South America Oceania S Africa Southeast Asia
## same 8 0 20 46
## distinct 0 12 0 20
This is perfectly uniform within the levels of each of the three factors, so random slopes are arguably not justified.
r/l distinction in L2 (r_l_distinction_L2):
table(webfield$r_l_distinction_L2, webfield$Language)
##
## AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PA PL PT RO
## same 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## distinct 9 20 18 20 30 98 16 43 28 34 18 18 46 15 42 26 50 28 21 8 52 42 30
##
## RU SE SR TH TR VA ZU
## same 0 0 0 0 0 0 0
## distinct 43 19 1 18 29 12 18
table(webfield$r_l_distinction_L2, webfield$Family)
##
## Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric IE
## same 0 0 0 0 0 2
## distinct 20 8 12 18 87 566
##
## Japanese Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
## same 1 0 0 0 0 0
## distinct 28 15 21 30 18 29
table(webfield$r_l_distinction_L2, webfield$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
## same 2 0 0 0 1
## distinct 511 82 104 20 49
##
## NE South America Oceania S Africa Southeast Asia
## same 0 0 0 0
## distinct 8 12 18 48
There is almost uniform, so random slopes are arguably not justified.
[r] is main allophone in L1 (trill_real_L1):
table(webfield$trill_real_L1, webfield$Language)
##
## AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT JP
## no 0 0 55 0 46 101 18 0 38 0 0 0 56 14 42 0 0 55
## yes 10 20 0 20 0 3 0 43 1 36 21 18 1 1 0 34 52 0
##
## KR PA PL PT RO RU SE SR TH TR VA ZU
## no 22 8 0 61 0 0 20 13 20 37 0 20
## yes 0 0 53 0 31 47 1 0 0 0 12 0
table(webfield$trill_real_L1, webfield$Family)
##
## Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric IE Japanese
## no 0 8 0 20 0 404 55
## yes 20 0 12 0 95 276 0
##
## Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
## no 14 22 46 20 37
## yes 1 0 0 0 0
table(webfield$trill_real_L1, webfield$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
## no 404 51 0 0 77
## yes 222 42 108 20 0
##
## NE South America Oceania S Africa Southeast Asia
## no 8 0 20 66
## yes 0 12 0 0
This is almost uniform within Language, but there is variation for the IE level of Family, and between 2 levels of Area (Europeand Greater Mesopotamia), but this is not enough to justify random slopes.
[r] is main allophone in L2 (trill_real_L2):
table(webfield$trill_real_L2, webfield$Language)
##
## AL AM BE BR CN DE DK EE EN ES FA FI FR GE GR HU IT JP KR PA PL PT RO RU
## no 5 1 9 0 29 48 13 5 16 22 16 8 28 3 28 18 30 26 21 8 34 19 21 32
## yes 4 19 10 20 1 50 3 38 13 12 2 10 18 12 14 8 20 3 0 0 18 23 9 11
##
## SE SR TH TR VA ZU
## no 11 1 18 28 0 16
## yes 8 0 0 1 12 2
table(webfield$trill_real_L2, webfield$Family)
##
## Afro-Asiatic Arawakan Austronesian Benue-Congo Finno-Ugric IE Japanese
## no 0 8 0 16 31 334 26
## yes 20 0 12 2 56 234 3
##
## Kartvelian Korean Sino-Tibetan Tai-Kadai Turkic
## no 3 21 29 18 28
## yes 12 0 1 0 1
table(webfield$trill_real_L2, webfield$Autotyp_Area)
##
## Europe Greater Mesopotamia Inner Asia N Africa N Coast Asia
## no 303 48 45 0 47
## yes 210 34 59 20 3
##
## NE South America Oceania S Africa Southeast Asia
## no 8 0 16 47
## yes 0 12 2 1
Here there is enough variation within the levels for all three factors, justifying the inclusion of random slopes.
Given these, we did the following:
if( !file.exists("./cached_results/b1_res_webfield.RData") )
{
# prior predictive checks:
# what priors we can set:
get_prior(Match ~ 1 + dataset +
(1 | Language) +
(1 | Family) +
(1 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'));
b_priors <- brms::brm(Match ~ 1 + dataset +
(1 | Language) +
(1 | Family) +
(1 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
sample_prior='only', # needed for prior predictive checks
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
grid.arrange(pp_check(b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is a bit extreme
# the model:
b1 <- brm(Match ~ 1 + dataset +
(1 | Language) +
(1 | Family) +
(1 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999999, max_treedepth=13));
summary(b1); mcmc_plot(b1, type="trace"); mcmc_plot(b1); # very decent
mcmc_plot(b1, type="trace");
bayestestR::hdi(b1, ci=0.95);
hypothesis(b1, "datasetfield = 0");
# posterior predictive checks
pp_check(b1, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values'),
pp_check(b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means'),
pp_check(b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values'),
ncol=1); # far from ideal but seems ok...
# save the model:
b1_res_webfield <- list("b_priors"=b_priors, "b1"=b1);
saveRDS(b1_res_webfield, "./cached_results/b1_res_webfield.RData", compress="xz");
} else
{
b1_res_webfield <- readRDS("./cached_results/b1_res_webfield.RData");
}
grid.arrange(pp_check(b1_res_webfield$b_priors, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Prior predictive distribution of minimum values'),
pp_check(b1_res_webfield$b_priors, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Prior predictive distribution of means'),
pp_check(b1_res_webfield$b_priors, type='stat', stat='max') + xlab('p(match)') + ggtitle('Prior predictive distribution of maximum values'),
ncol=1); # seems fine but our value is a bit extreme
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 29. Prior predictive checks. Figure generated
using R version
4.3.3 (2024-02-29)
and it can be seen that they are quite ok, even if the observed p(match) is a bit extreme, but still within what the prior distribution covers.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b1_res_webfield$b1);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + dataset + (1 | Language) + (1 | Family) + (1 | Autotyp_Area)
## Data: webfield (Number of observations: 1030)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.40 0.36 0.02 1.31 1.00 6971 8935
##
## ~Family (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.34 0.29 0.01 1.05 1.00 7936 9857
##
## ~Language (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.57 0.19 0.23 0.98 1.00 6446 7944
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.01 0.33 1.39 2.71 1.00 10336 10072
## datasetfield 1.64 0.68 0.39 3.06 1.00 10442 10343
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b1_res_webfield$b1, type="trace");
## No divergences to plot.
Figure 30. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks seems to be acceptable:
pp_check(b1_res_webfield$b1, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 31. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(pp_check(b1_res_webfield$b1, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b1_res_webfield$b1, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b1_res_webfield$b1, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 32. Posterior predictive distribution. Figure
generated using R
version 4.3.3 (2024-02-29)
Fixed effects:
mcmc_plot(b1_res_webfield$b1);
Figure 33. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
plot_model(b1_res_webfield$b1, type="emm", terms=c("dataset"));
Figure 34. Fixed effects. Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 2.01 95%CI [1.39, 2.71] on the log odds ratio (LOR) scale, which is clearly » 0; formal hypothesis testing against 0 is p(α>0)=1.00*; this translates into a probability of a match p(match)=88.2% 95%CI [80.0%, 93.8%] » 50%. (N.B. this is for dataset == “web”; please below for details.)
Dataset (experiment) The slope of dataset is 1.64 95%CI [0.39, 3.06] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(β>0)=0.99*. In terms of probabilities, the field experiment (dataset == “field”) results in an increase by 9.3% in the probability of a match from the 88.2% [80.0%, 93.8%] in the web experiment, to 97.5% [91.5%, 99.4%] in the field experiment.
Random effects:
#get_variables(b1_res_webfield$b1);
grid.arrange(
# 1 | Language:
b1_res_webfield$b1 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Family:
b1_res_webfield$b1 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Autotyp_Area:
b1_res_webfield$b1 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
nrow=1);
Figure 35. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match, while controlling for Language, Family and Area, is much higher than the conservative change level of 50%, being estimated at more than 88% [80%, 94%] in the web experiment, and more than 97% [92%, 99%] in the field experiment, suggesting that there is a very strong tendency across our participants to associate [r] with the wavy line and [l] with the straight line.
if( !file.exists("./cached_results/b2_res_webfield.RData") )
{
# the model:
b2 <- brm(Match ~ 1 + dataset + # dataset
r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + # main effects
dataset:r_l_distinction_L1 + dataset:r_l_distinction_L2 + dataset:trill_real_L1 + dataset:trill_real_L2 + # potential interactions with dataset
(1 + trill_real_L2 | Language) +
(1 + trill_real_L2 | Family) +
(1 + trill_real_L2 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"),
brms::set_prior("lkj(2)", class="cor")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b2); mcmc_plot(b2, type="trace"); mcmc_plot(b2); # very decent
mcmc_plot(b2, type="trace");
bayestestR::hdi(b2, ci=0.95);
hypothesis(b2, c("datasetfield = 0",
"r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0",
"datasetfield:r_l_distinction_L1distinct = 0", "datasetfield:r_l_distinction_L2distinct = 0", "datasetfield:trill_real_L1yes = 0", "datasetfield:trill_real_L2yes = 0"));
# posterior predictive checks
pp_check(b2, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b2, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1); # seems perfect
# the interactions of dataset do not seem to matter at all, so let's remove them:
b3 <- brm(Match ~ 1 + dataset + # dataset
r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + # main effects
(1 + trill_real_L2 | Language) +
(1 + trill_real_L2 | Family) +
(1 + trill_real_L2 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b"),
brms::set_prior("lkj(2)", class="cor")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b3); mcmc_plot(b3, type="trace"); mcmc_plot(b3); # very decent
mcmc_plot(b3, type="trace");
bayestestR::hdi(b3, ci=0.95);
hypothesis(b3, c("datasetfield = 0", "r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0"));
# posterior predictive checks
pp_check(b3, ndraws=100) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
grid.arrange(pp_check(b3, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b3, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b3, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1); # seems perfect
# save the model:
b2_res_webfield <- list("b2"=b2, "b3"=b3);
saveRDS(b2_res_webfield, "./cached_results/b2_res_webfield.RData", compress="xz");
} else
{
b2_res_webfield <- readRDS("./cached_results/b2_res_webfield.RData");
}
We fist included the interactions between dataset and the [r]-related predictors, but they are not relevant (it is interesting to note that even in this model, trill_real_L1 has a significant effect):
bayestestR::hdi(b2_res_webfield$b2, ci=0.95);
hypothesis(b2_res_webfield$b2, c("datasetfield = 0",
"r_l_distinction_L1distinct = 0", "r_l_distinction_L2distinct = 0", "trill_real_L1yes = 0", "trill_real_L2yes = 0",
"datasetfield:r_l_distinction_L1distinct = 0", "datasetfield:r_l_distinction_L2distinct = 0", "datasetfield:trill_real_L1yes = 0", "datasetfield:trill_real_L2yes = 0"));
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (datasetfield) = 0 2.49 2.95 -2.31 9.31 0.80
## 2 (r_l_distinction_... = 0 0.49 0.93 -1.47 2.25 2.55
## 3 (r_l_distinction_... = 0 -1.29 2.59 -7.11 2.59 1.26
## 4 (trill_real_L1yes) = 0 -0.99 0.40 -1.81 -0.24 0.27
## 5 (trill_real_L2yes) = 0 0.38 0.93 -1.27 2.47 3.52
## 6 (datasetfield:r_l... = 0 1.50 2.73 -3.17 7.78 1.05
## 7 (datasetfield:r_l... = 0 2.24 2.86 -2.46 8.74 0.88
## 8 (datasetfield:tri... = 0 0.80 2.78 -4.31 6.84 1.09
## 9 (datasetfield:tri... = 0 0.95 2.86 -4.09 6.88 1.12
## Post.Prob Star
## 1 0.44
## 2 0.72
## 3 0.56
## 4 0.21 *
## 5 0.78
## 6 0.51
## 7 0.47
## 8 0.52
## 9 0.53
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
so we removed them, keeping only the main effects of the predictors (including dataset).
Priors are as before.
Convergence The model converges well (high ESS and \(\hat{R} = 1.0\)):
summary(b2_res_webfield$b3);
## Family: bernoulli
## Links: mu = logit
## Formula: Match ~ 1 + dataset + r_l_distinction_L1 + r_l_distinction_L2 + trill_real_L1 + trill_real_L2 + (1 + trill_real_L2 | Language) + (1 + trill_real_L2 | Family) + (1 + trill_real_L2 | Autotyp_Area)
## Data: webfield (Number of observations: 855)
## Draws: 4 chains, each with iter = 10000; warmup = 4000; thin = 2;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~Autotyp_Area (Number of levels: 9)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.57 0.50 0.02 1.89 1.00
## sd(trill_real_L2yes) 1.07 0.93 0.04 3.48 1.00
## cor(Intercept,trill_real_L2yes) -0.03 0.45 -0.83 0.81 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 6674 9320
## sd(trill_real_L2yes) 7321 8492
## cor(Intercept,trill_real_L2yes) 9759 9994
##
## ~Family (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.93 0.56 0.06 2.19 1.00
## sd(trill_real_L2yes) 0.95 0.89 0.03 3.22 1.00
## cor(Intercept,trill_real_L2yes) -0.16 0.45 -0.88 0.74 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5023 6050
## sd(trill_real_L2yes) 7748 9334
## cor(Intercept,trill_real_L2yes) 10006 10648
##
## ~Language (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.38 0.27 0.02 0.99 1.00
## sd(trill_real_L2yes) 0.87 0.48 0.08 1.94 1.00
## cor(Intercept,trill_real_L2yes) -0.10 0.44 -0.84 0.75 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 5731 7352
## sd(trill_real_L2yes) 5253 6479
## cor(Intercept,trill_real_L2yes) 6527 9138
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 3.27 2.57 -0.82 9.14 1.00 8579
## datasetfield 4.22 2.58 0.92 10.16 1.00 8095
## r_l_distinction_L1distinct 0.54 0.93 -1.45 2.33 1.00 8002
## r_l_distinction_L2distinct -1.21 2.47 -6.85 2.59 1.00 8657
## trill_real_L1yes -0.98 0.41 -1.81 -0.23 1.00 8916
## trill_real_L2yes 0.47 0.95 -1.12 2.70 1.00 8081
## Tail_ESS
## Intercept 8945
## datasetfield 6470
## r_l_distinction_L1distinct 8899
## r_l_distinction_L2distinct 8438
## trill_real_L1yes 9618
## trill_real_L2yes 8436
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_plot(b2_res_webfield$b3, type="trace");
## No divergences to plot.
Figure 36. Trace of the fitting process. Figure
generated using R
version 4.3.3 (2024-02-29)
Posterior predictive checks seems to be perfectly fine:
pp_check(b2_res_webfield$b3, ndraws=500) + xlab('p(match)') + ggtitle('Posterior predictive density overlay');
Figure 37. Posterior predictive density overlay for
n=500 draws. Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(pp_check(b2_res_webfield$b3, type='stat', sta ='min') + xlab('p(match)') + ggtitle('Posterior predictive distribution of minimum values') + xlim(0,1),
pp_check(b2_res_webfield$b3, type='stat', stat='mean') + xlab('p(match)') + ggtitle('Posterior predictive distribution of means') + xlim(0,1),
pp_check(b2_res_webfield$b3, type='stat', stat='max') + xlab('p(match)') + ggtitle('Posterior predictive distribution of maximum values') + xlim(0,1),
ncol=1);
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## Using all posterior draws for ppc type 'stat' by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Figure 38. Posterior predictive distribution (the
x-axis scale was forced to [0,1]). Figure generated using R version 4.3.3
(2024-02-29)
Fixed effects:
mcmc_plot(b2_res_webfield$b3);
Figure 39. Fixed effects’ posterior distributions.
Figure generated using R version 4.3.3
(2024-02-29)
grid.arrange(plot_model(b2_res_webfield$b3, type="emm", terms=c("dataset")),
plot_model(b2_res_webfield$b3, type="emm", terms=c("r_l_distinction_L1")),
plot_model(b2_res_webfield$b3, type="emm", terms=c("r_l_distinction_L2")),
plot_model(b2_res_webfield$b3, type="emm", terms=c("trill_real_L1")),
plot_model(b2_res_webfield$b3, type="emm", terms=c("trill_real_L2")),
ncol=2);
Figure 40. Fixed effects. Figure generated using R version 4.3.3
(2024-02-29)
Intercept The intercept is 3.27 95%CI [-0.82, 9.14] on the log odds ratio (LOR) scale, which is > 0 but has very large uncertainties so that the 95%CI includes 0; formal hypothesis testing against 0 is p(α>0)=0.93; this translates into a probability of a match p(match)=96.4% 95%CI [30.6%, 100.0%] which includes 50% in the 95%CI.
Dataset (experiment) The slope of dataset is 4.22 95%CI [0.92, 10.16] on the log odds ratio (LOR) scale, which is clearly > 0; formal hypothesis testing against 0 is p(β>0)=1.00*. In terms of probabilities, the field experiment (dataset == “field”) results in an increase by 3.6% in the probability of a match from the 96.4% [30.6%, 100.0%] in the web experiment, to 99.9% [87.8%, 100.0%] in the field experiment.
Are [r] and [l] distinct in the L1(s) spoken by the participant? The slope of r_l_distinction_L1 is 0.54 95%CI [-1.45, 2.33] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β>0)=0.75.
Are [r] and [l] distinct in the L2(s) spoken by the participant? The slope of r_l_distinction_L2 is -1.21 95%CI [-6.85, 2.59] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 is p(β<0)=0.67.
Is [r] the main allophone in the L1(s) spoken by the participant? The slope of trill_real_L1 is -0.98 95%CI [-1.81, -0.23] on the log odds ratio (LOR) scale, which is clearly < 0; formal hypothesis testing against 0 is p(β<0)=0.99*. In terms of probabilities, having [r] as the main allophone (trill_real_L1 == “yes”) results in a decrease by 5.6% [0.003%, 23.9%] in the probability of a match from the 96.4% [30.6%, 100.0%] when [r] is not the primary allophone to 90.8% [13.7%, 100.0%] when it is.
Is [r] the main allophone in the L2(s) spoken by the participant? The slope of trill_real_L2 is 0.47 95%CI [-1.12, 2.70] on the log odds ratio (LOR) scale, which includes 0; formal hypothesis testing against 0 with p(β>0)=0.68.
Random effects:
#get_variables(b2_res_webfield$b3);
grid.arrange(
# 1 | Language:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Language[Language, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Language") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Family:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Family[Family, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Family") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# 1 | Autotyp_Area:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, Order]) %>%
filter(Order == "Intercept") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("1 | Area") + xlab("Intercept") + ylab(NULL) + xlim(-2,2) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Language:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Language[Language, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Language) %>%
ggplot(aes(y = Language, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Language") + xlab("Slope") + ylab(NULL) + xlim(-4, 4) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Family:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Family[Family, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Family) %>%
ggplot(aes(y = Family, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Family") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
# trill_real_L2 | Autotyp_Area:
b2_res_webfield$b3 %>%
spread_draws(b_Intercept, r_Autotyp_Area[Autotyp_Area, trill_real_L2]) %>%
filter(trill_real_L2 == "trill_real_L2yes") %>%
mutate(condition_mean = r_Autotyp_Area) %>%
ggplot(aes(y = Autotyp_Area, x = condition_mean, fill = after_stat(x < 0))) +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
geom_vline(xintercept = 0, linetype = "dashed") +
ggtitle("trill_real_L2 | Area") + xlab("Slope") + ylab(NULL) + xlim(-3, 3) +
theme(legend.position="bottom") +
NULL,
nrow=2);
Figure 41. Posterior estimaes of the random effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
Interpretation: The overall probability of a “perfect” match in the online experiment when adding all the relevant [r]-related predictors for the L1(s) and L2(s) known to the participants, while controlling for Language, Family and Area, seem to still be higher than the conservative change level of 50%, being estimated as about 96% but with a much wider 95% credible interval of about [31%, 100.0%] which now includes 50%. If [r] is the main allophone in the L1(s) known by the participant or not clearly matters, in that those participants that have an L1 with [r] as the main allophone have a slightly smaller probability of a perfect match of about 91% [14%, 100%] than the others. Also, there is a clear difference between the two datasets, with the web participants having a lower probability of a match than the field participants by about 4%, to 100% [88%, 100%].
Here we build the various plots for the paper.
p <- gather_draws(b2_res_web$b2, `b_.*`, regex = TRUE) %>%
# remove the b_ for plotting
#mutate(.variable = gsub("b_","", .variable)) %>%
ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
xlim(-10, 10) +
scale_y_discrete(labels=c("b_Intercept" = "Intercept α",
"b_Order1" = "order ([r] first vs [l] first)",
"b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
"b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
"b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
"b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?")) +
ylab(NULL) + xlab("Estimate") +
theme_timo + # Tweak cosmetics:
theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
r = 0, l = 0),
face = 'bold', size = 14),
axis.text.x = element_text(face = 'bold', size = 12),
axis.text.y = element_text(face = 'bold', size = 12))+
NULL;
# Show and save:
p;
Figure 42. Posterior estimaes of the fixed effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/figure_model_fixedeffects_web.tif", width = 8, height = 4, compression="lzw", dpi=600);
# take2: generate an exhaustive new dataset:
# starting from the actual language coding:
languages_data <- read.csv("./data/languages_data.csv", sep=","); # read the original language info
newdata_web <- data.frame(Name = tolower(unique(web$Name)));
newdata_web <- merge(newdata_web, languages_data, by.x="Name", by.y="Languages", all.x=TRUE, all.y=FALSE);
# L1 coding:
names(newdata_web)[3:5] <- c("r_l_distinction_L1", "trill_real_L1", "trill_occ_L1");
newdata_web$r_l_distinction_L1 <- factor(c("same", "distinct")[newdata_web$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
newdata_web$trill_real_L1 <- factor(c("no", "yes")[newdata_web$trill_real_L1 + 1], levels=c("no", "yes"));
newdata_web$trill_occ_L1 <- factor(c("no", "yes")[newdata_web$trill_occ_L1 + 1], levels=c("no", "yes"));
# family and area:
newdata_web <- merge(newdata_web, unique(web[,c("Language", "Name", "Family", "Autotyp_Area")]) %>% mutate(Name = tolower(Name)), by="Name", all.x=TRUE, all.y=FALSE);
# add it to the both newdata:
newdata_webfield <- newdata_web[ ,c("Name", "Language", "Family", "Autotyp_Area", "r_l_distinction_L1", "trill_real_L1", "trill_occ_L1")];
# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_web.RData") )
{
b_plotting_web <- brm(Match ~ 1 +
r_l_distinction_L1 + trill_real_L1 +
(1 | Language) +
(1 | Family) +
(1 | Autotyp_Area),
data=web,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b_plotting_web); mcmc_plot(b_plotting_web, type="trace"); mcmc_plot(b_plotting_web); # very decent
bayestestR::hdi(b_plotting_web, ci=0.95);
# save the model:
saveRDS(b_plotting_web, "./cached_results/b_plotting_web.RData", compress="xz");
} else
{
b_plotting_web <- readRDS("./cached_results/b_plotting_web.RData");
}
fit <- fitted(b_plotting_web, newdata=newdata_web, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_web <- cbind(newdata_web, fit);
# Order predictions by descriptive average:
newdata_web <- arrange(newdata_web, fit);
newdata_web <- mutate(newdata_web, Language = factor(Language, levels = newdata_web$Language));
# How many languages are over 0.5? (will be reported in paper)
sum(newdata_web$lwr > 0.5);
# Finally, add the averages to the plot:
newdata_web$avg <- web_avg[match(newdata_web$Language, web_avg$Language), ]$M;
# Match language names into there and order:
newdata_web$Language <- web[match(newdata_web$Language, web$Language), ]$Name;
newdata_web[newdata_web$Language == 'Chinese', ]$Language <- 'Mandarin Chinese';
newdata_web <- mutate(newdata_web, Language = factor(Language, levels = newdata_web$Language));
# Setup the plot:
# Aesthetics and geom:
plot_web <- newdata_web %>%
ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
#scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
geom_point(size = 5) +
geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') +
geom_point(aes(y = avg), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.75) +
ylim(0.5,1.00) +
labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
#ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) +
theme_timo +
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
axis.text.y = element_text(face = 'bold', size = 10),
axis.title = element_text(face = 'bold', size = 12),
axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
legend.text = element_text(size = 10),
legend.title = element_blank(),
legend.position = c(0.95, 0.10),
legend.justification = c('right', 'bottom'))
# Save:
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.pdf", width = 12, height = 6);
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.jpg", width = 12, height = 6);
ggsave(plot = plot_web, filename = "./plots/figure_model_languages_web.tif", width = 12, height = 6, compression="lzw", dpi=600);
plot_web;
Figure 43. The distribution of the proportion of
matches per language. Please note that this was obtained by
fitting a model that includes only the relevant L1 characteristics,
Match ~ 1 + r_l_distinction_L1 + trill_real_L1 + (1 | Language) + (1 | Family) + (1 | Autotyp_Area).
Figure generated using R version 4.3.3
(2024-02-29)
p <- gather_draws(b2_res_field$b2, `b_.*`, regex = TRUE) %>%
# remove the b_ for plotting
#mutate(.variable = gsub("b_","", .variable)) %>%
ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
xlim(-10, 20) +
scale_y_discrete(labels=c("b_Intercept" = "Intercept α",
"b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
"b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
"b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
"b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?",
"b_trill_occ_L1yes" = "is [r] an allophone in L1(s)?")) +
ylab(NULL) + xlab("Estimate") +
theme_timo + # Tweak cosmetics:
theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
r = 0, l = 0),
face = 'bold', size = 14),
axis.text.x = element_text(face = 'bold', size = 12),
axis.text.y = element_text(face = 'bold', size = 12))+
NULL;
# Show and save:
p;
Figure 44. Posterior estimaes of the fixed effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
ggsave(plot = p, filename = "./plots/full_model_field_coefficients.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_field_coefficients.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_field_coefficients.tif", width = 8, height = 4, compression="lzw", dpi=600);
# take2: generate an exhaustive new dataset:
# starting from the actual language coding:
languages_data <- read.csv("./data/languages_data.csv", sep=","); # read the original language info
newdata_field <- data.frame(Name = tolower(unique(field$Name)));
newdata_field$Name2 <- newdata_field$Name;
newdata_field$Name[ newdata_field$Name == "brazilian portuguese" ] <- "portuguese";
newdata_field$Name[ newdata_field$Name == "english uk" ] <- "english";
newdata_field$Name[ newdata_field$Name == "tashlhiyt berber" ] <- "berber";
newdata_field <- merge(newdata_field, languages_data, by.x="Name", by.y="Languages", all.x=TRUE, all.y=FALSE);
# manually fix some entries:
# L1 coding:
names(newdata_field)[4:6] <- c("r_l_distinction_L1", "trill_real_L1", "trill_occ_L1");
newdata_field$r_l_distinction_L1 <- factor(c("same", "distinct")[newdata_field$r_l_distinction_L1 + 1], levels=c("same", "distinct"));
newdata_field$trill_real_L1 <- factor(c("no", "yes")[newdata_field$trill_real_L1 + 1], levels=c("no", "yes"));
newdata_field$trill_occ_L1 <- factor(c("no", "yes")[newdata_field$trill_occ_L1 + 1], levels=c("no", "yes"));
# family and area:
newdata_field <- merge(newdata_field, unique(field[,c("Language", "Name", "Family", "Autotyp_Area")]) %>% mutate(Name = tolower(Name)), by.x="Name2", by.y="Name", all.x=TRUE, all.y=FALSE);
# add it to the both newdata:
newdata_webfield <- rbind(newdata_webfield,
newdata_field[,c("Name", "Language", "Family", "Autotyp_Area", "r_l_distinction_L1", "trill_real_L1", "trill_occ_L1")]);
# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_field.RData") )
{
b_plotting_field <- brm(Match ~ 1 +
r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 +
(1 | Language),
data=field,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b_plotting_field); mcmc_plot(b_plotting_field, type="trace"); mcmc_plot(b_plotting_field); # very decent
bayestestR::hdi(b_plotting_field, ci=0.95);
# save the model:
saveRDS(b_plotting_field, "./cached_results/b_plotting_field.RData", compress="xz");
} else
{
b_plotting_field <- readRDS("./cached_results/b_plotting_field.RData");
}
fit <- fitted(b_plotting_field, newdata=newdata_field, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_field <- cbind(newdata_field, fit);
# Order predictions by descriptive average:
newdata_field <- arrange(newdata_field, fit);
newdata_field <- mutate(newdata_field, Language = factor(Language, levels = newdata_field$Language));
# How many languages are over 0.5? (will be reported in paper)
sum(newdata_field$lwr > 0.5);
# Finally, add the averages to the plot:
newdata_field$avg <- field_avg[match(newdata_field$Language, field_avg$Language), ]$M;
# Match language names into there and order:
newdata_field$Language <- field[match(newdata_field$Language, field$Language), ]$Name;
newdata_field <- mutate(newdata_field, Language = factor(Language, levels = newdata_field$Language));
# Setup the plot:
# Aesthetics and geom:
plot_field <- newdata_field %>%
ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
#scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
geom_point(size = 5) +
geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') +
geom_point(aes(y = avg), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.75) +
ylim(0.5,1.00) +
labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
#ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) +
theme_timo +
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
axis.text.y = element_text(face = 'bold', size = 10),
axis.title = element_text(face = 'bold', size = 12),
axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
legend.text = element_text(size = 10),
legend.title = element_blank(),
legend.position = c(0.95, 0.10),
legend.justification = c('right', 'bottom'))
# Save:
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.pdf", width = 6, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.jpg", width = 6, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.tif", width = 6, height = 6, compression="lzw", dpi=600);
plot_field;
Figure 45. The distribution of the proportion of
matches per language. Please note that this was obtained by
fitting a model that includes only the relevant L1 characteristics,
Match ~ 1 + r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 + (1 | Language).
Figure generated using R version 4.3.3
(2024-02-29)
p <- gather_draws(b2_res_webfield$b3, `b_.*`, regex = TRUE) %>%
# remove the b_ for plotting
#mutate(.variable = gsub("b_","", .variable)) %>%
ggplot(aes(x = .value, y = .variable, fill = after_stat(x < 0))) +
geom_vline(xintercept=0.0, linetype="dotted", color="gray30") +
stat_halfeye(alpha=0.75, point_interval=median_qi, .width=c(.50, .95)) +
scale_fill_manual(values = c("lightsalmon", "skyblue"), name=NULL, labels=c("≥ 0", "< 0")) +
xlim(-10, 10) +
scale_y_discrete(labels=c("b_Intercept" = "Intercept α",
"b_datasetfield" = "the field vs the web dataset",
"b_r_l_distinction_L1distinct" = "are [r] and [l] contrasting in L1(s)?",
"b_trill_real_L1yes" = "is [r] the main allophone in L1(s)?",
"b_r_l_distinction_L2distinct" = "are [r] and [l] contrasting in L2(s)?",
"b_trill_real_L2yes" = "is [r] the main allophone in L2(s)?")) +
ylab(NULL) + xlab("Estimate") +
theme_timo + # Tweak cosmetics:
theme(axis.title.x = element_text(margin = margin(t = 12, b = 0,
r = 0, l = 0),
face = 'bold', size = 14),
axis.text.x = element_text(face = 'bold', size = 12),
axis.text.y = element_text(face = 'bold', size = 12))+
NULL;
# Show and save:
p;
Figure 46. Posterior estimaes of the fixed effects
showing the median (dot), 50% (thick line) and 95% (thin lines)
quantiles. Also showing 0 as the vertical dotted line and negative
(blue) vs positive (red) values. Figure generated using R version 4.3.3
(2024-02-29)
ggsave(plot = p, filename = "./plots/full_model_both_coefficients.pdf", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_both_coefficients.jpg", width = 8, height = 4);
ggsave(plot = p, filename = "./plots/full_model_both_coefficients.tif", width = 8, height = 4, compression="lzw", dpi=600);
# take2: generate an exhaustive new dataset:
newdata_webfield$dataset <- factor(c(rep("web", nrow(newdata_web)), rep("field", nrow(newdata_field))), levels=c("web", "field"));
# fit a new model considering only the L1 info:
if( !file.exists("./cached_results/b_plotting_webfield.RData") )
{
b_plotting_webfield <- brm(Match ~ 1 + dataset +
r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 +
(1 | Language) +
(1 | Family) +
(1 | Autotyp_Area),
data=webfield,
family=bernoulli(link='logit'),
prior=c(brms::set_prior("student_t(5, 0, 2.5)", class="b")),
save_pars=save_pars(all=TRUE), # needed for Bayes factors
sample_prior=TRUE, # needed for hypotheses tests
seed=998, cores=brms_ncores, iter=10000, warmup=4000, thin=2, control=list(adapt_delta=0.999, max_treedepth=13));
summary(b_plotting_webfield); mcmc_plot(b_plotting_webfield, type="trace"); mcmc_plot(b_plotting_webfield); # very decent
bayestestR::hdi(b_plotting_webfield, ci=0.95);
# save the model:
saveRDS(b_plotting_webfield, "./cached_results/b_plotting_webfield.RData", compress="xz");
} else
{
b_plotting_webfield <- readRDS("./cached_results/b_plotting_webfield.RData");
}
fit <- fitted(b_plotting_webfield, newdata=newdata_webfield, re_formula=NULL, robust=TRUE);
colnames(fit) = c('fit', 'se', 'lwr', 'upr');
newdata_webfield <- cbind(newdata_webfield, fit);
# Order predictions by descriptive average:
newdata_webfield <- arrange(newdata_webfield, fit);
#newdata_webfield <- mutate(newdata_webfield, Language = factor(Language, levels = newdata_webfield$Language));
# How many languages are over 0.5? (will be reported in paper)
sum(newdata_webfield$lwr > 0.5);
# Finally, add the averages to the plot:
newdata_webfield$avg[ newdata_webfield$dataset == "web" ] <- web_avg[match(newdata_webfield$Language[ newdata_webfield$dataset == "web" ], web_avg$Language), ]$M;
newdata_webfield$avg[ newdata_webfield$dataset == "field" ] <- field_avg[match(newdata_webfield$Language[ newdata_webfield$dataset == "field" ], field_avg$Language), ]$M;
# Match language names into there and order:
newdata_webfield$Language[ newdata_webfield$dataset == "web" ] <- paste0(web[match(newdata_webfield$Language[ newdata_webfield$dataset == "web" ], web$Language), ]$Name, " (W)");
newdata_webfield$Language[ newdata_webfield$dataset == "field" ] <- paste0(field[match(newdata_webfield$Language[ newdata_webfield$dataset == "field" ], field$Language), ]$Name, " (F)");
newdata_webfield <- mutate(newdata_webfield, Language = factor(Language, levels = newdata_webfield$Language));
# Setup the plot:
# Aesthetics and geom:
plot_field <- newdata_webfield %>%
ggplot(aes(x = Language, col = r_l_distinction_L1, y = fit, ymin = lwr, ymax = upr, shape = trill_real_L1)) +
#scale_shape_manual(values = c("no" = 1, "yes" = 3)) +
geom_errorbar(aes(col = r_l_distinction_L1), linewidth = 1.0, width = 0.3) +
geom_point(size = 5) +
geom_hline(yintercept = 0.5, linetype = 2, linewidth = 1.5, col = 'grey') +
geom_point(aes(y = avg, fill=dataset), col = 'black', shape = 23, size = 5, stroke = 1.0, alpha = 0.25) +
ylim(0.5,1.00) +
labs(x = '', y = 'Proportion of congruent responses') + # Axis labels:
#ggtitle('Posterior medians and descriptive averages of congruent responses\nby language and R/L contrast (color) and [r] as primary rhotic in L1 (shape)') +
scale_color_manual(values = c("same"=colorBlindBlack8[2], "distinct"=colorBlindBlack8[3]), labels=c("same"="no [r]/[l] contrast", "distinct"="[r]/[l] contrast")) + # Tweak cosmetics:
scale_shape(labels=c("yes"="[r] is main allophone", "no"="[r] is not main allophone")) +
scale_fill_manual(values = c("web"="white", "field"="black"), labels=c("web"="web experiment", "field"="field experiment")) +
theme_timo +
theme(axis.text.x = element_text(angle = 45, hjust = 1, face = 'bold', size = 10),
axis.text.y = element_text(face = 'bold', size = 10),
axis.title = element_text(face = 'bold', size = 12),
axis.title.y = element_text(face = 'bold', size = 12, margin = margin(t = 0, r = 35, b = 0, l = 0)),
plot.title = element_text(face = 'bold', size = 14, margin = margin(t = 0, r = 0, b = 30, l = 0)),
legend.text = element_text(size = 10),
legend.title = element_blank(),
legend.position = c(0.95, 0.10),
legend.justification = c('right', 'bottom'))
# Save:
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.pdf", width = 18, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.jpg", width = 18, height = 6);
ggsave(plot = plot_field, filename = "./plots/figure_model_languages_field.tif", width = 18, height = 6, compression="lzw", dpi=600);
plot_field;
Figure 47. The distribution of the proportion of
matches per language. Please note that this was obtained by
fitting a model that includes only the relevant L1 characteristics,
Match ~ 1 + dataset + r_l_distinction_L1 + trill_real_L1 + trill_occ_L1 + (1 | Language) + (1 | Family) + (1 | Autotyp_Area).
Figure generated using R version 4.3.3
(2024-02-29)
It is clear that across the two datasets our participants have a strong tendency to associate [r] with the zigzagy line and [l] with the straight line well above the conservative 50% chance level. Presenting [r] first seems to increase the probability of a match in the online experiment, compared to presenting [l] first. The field experiment has an overall higher probability of a match than the online experiment. In the online experiment and in the combined datasets, having [r] as the primary allophone in the the L1(s) spoken by a participant slightly decreases the probability of a match.
if( require(benchmarkme) )
{
# CPU:
cpu_info <- benchmarkme::get_cpu();
if( is.null(cpu_info) || all(is.na(cpu_info)) )
{
cat("**CPU:** unknown.\n\n");
} else
{
if( !is.null(cpu_info$model_name) && !is.na(cpu_info$model_name) )
{
cat(paste0("**CPU:** ",cpu_info$model_name));
if( !is.null(cpu_info$no_of_cores) && !is.na(cpu_info$no_of_cores) )
{
cat(paste0(" (",cpu_info$no_of_cores," threads)"));
}
cat("\n\n");
} else
{
cat("**CPU:** unknown.\n\n");
}
}
# RAM:
ram_info <- benchmarkme::get_ram();
if( is.null(ram_info) || is.na(ram_info) )
{
cat("**RAM (memory):** unknown.\n\n");
} else
{
cat("**RAM (memory):** "); print(ram_info); cat("\n");
}
} else
{
cat("**RAM (memory):** cannot get info (try installing package 'benchmarkme').\n\n");
}
CPU: Apple M3 (8 threads)
RAM (memory): 17.2 GB
pander::pander(sessionInfo());
R version 4.3.3 (2024-02-29)
Platform: aarch64-apple-darwin20 (64-bit)
locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8
attached base packages: parallel, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: benchmarkme(v.1.0.8), knitr(v.1.48), sjPlot(v.2.8.16), ggpubr(v.0.6.0), gridExtra(v.2.3), bayestestR(v.0.14.0), performance(v.0.12.2), lme4(v.1.1-35.5), Matrix(v.1.6-5), tidybayes(v.3.0.6), ggdist(v.3.3.2), brms(v.2.21.0), Rcpp(v.1.0.13), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)
loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), tensorA(v.0.36.2.1), rstudioapi(v.0.16.0), jsonlite(v.1.8.8), datawizard(v.0.10.0), magrittr(v.2.0.3), TH.data(v.1.1-2), estimability(v.1.5.1), farver(v.2.1.2), nloptr(v.2.1.1), rmarkdown(v.2.27), ragg(v.1.3.2), vctrs(v.0.6.5), minqa(v.1.2.7), rstatix(v.0.7.2), htmltools(v.0.5.8.1), distributional(v.0.4.0), curl(v.5.2.1), haven(v.2.5.4), broom(v.1.0.6), sjmisc(v.2.8.10), sass(v.0.4.9), StanHeaders(v.2.32.9), bslib(v.0.8.0), plyr(v.1.8.9), sandwich(v.3.1-0), emmeans(v.1.10.2), zoo(v.1.8-12), cachem(v.1.1.0), TMB(v.1.9.14), iterators(v.1.0.14), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), sjlabelled(v.1.2.0), R6(v.2.5.1), fastmap(v.1.2.0), digest(v.0.6.36), numDeriv(v.2016.8-1.1), colorspace(v.2.1-1), textshaping(v.0.4.0), labeling(v.0.4.3), fansi(v.1.0.6), timechange(v.0.3.0), httr(v.1.4.7), abind(v.1.4-5), mgcv(v.1.9-1), compiler(v.4.3.3), pander(v.0.6.5), doParallel(v.1.0.17), bit64(v.4.0.5), withr(v.3.0.1), backports(v.1.5.0), inline(v.0.3.19), carData(v.3.0-5), QuickJSR(v.1.3.1), pkgbuild(v.1.4.4), highr(v.0.11), ggsignif(v.0.6.4), MASS(v.7.3-60.0.1), sjstats(v.0.19.0), loo(v.2.7.0), tools(v.4.3.3), glue(v.1.7.0), nlme(v.3.1-165), grid(v.4.3.3), checkmate(v.2.3.2), reshape2(v.1.4.4), generics(v.0.1.3), glmmTMB(v.1.1.9), gtable(v.0.3.5), tzdb(v.0.4.0), hms(v.1.1.3), car(v.3.1-2), utf8(v.1.2.4), foreach(v.1.5.2), pillar(v.1.9.0), vroom(v.1.6.5), posterior(v.1.5.0), benchmarkmeData(v.1.0.4), splines(v.4.3.3), lattice(v.0.22-6), survival(v.3.7-0), bit(v.4.0.5), tidyselect(v.1.2.1), arrayhelpers(v.1.1-0), V8(v.4.4.2), stats4(v.4.3.3), xfun(v.0.46), bridgesampling(v.1.1-2), matrixStats(v.1.3.0), rstan(v.2.32.6), stringi(v.1.8.4), yaml(v.2.3.10), boot(v.1.3-30), evaluate(v.0.24.0), codetools(v.0.2-20), cli(v.3.6.3), RcppParallel(v.5.1.8), systemfonts(v.1.1.0), xtable(v.1.8-4), munsell(v.0.5.1), jquerylib(v.0.1.4), ggeffects(v.1.7.0), coda(v.0.19-4.1), svUnit(v.1.0.6), rstantools(v.2.4.0), bayesplot(v.1.11.1), Brobdingnag(v.1.2-9), viridisLite(v.0.4.2), mvtnorm(v.1.2-5), scales(v.1.3.0), insight(v.0.20.2), crayon(v.1.5.3), rlang(v.1.1.4) and multcomp(v.1.4-26)